# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/7/24

library(ggrepel)
library(plyr)
library(gridExtra)
library(scales)
library(egg)
library(optparse)
library(grid)
library(cowplot)
library(tidyverse)
library(ggpubr)
library(FSA)
library(magrittr)

createWhenNoExist <- function(f) {
  !dir.exists(f) && dir.create(f)
}

asplit <- function(x, MARGIN) {
  dl <- length(dim(x))
  if (!dl) stop("dim(x) must have a positive length")
  if (is.object(x)) x <- if (dl == 2) as.matrix(x) else as.array(x)
  d <- dim(x)
  dn <- dimnames(x)
  ds <- seq_len(dl)
  if (is.character(MARGIN)) {
    if (is.null(dnn <- names(dn))) stop("\'x\' must have named dimnames")
    MARGIN <- match(MARGIN, dnn)
    if (anyNA(MARGIN)) stop("not all elements of \'MARGIN\' are names of dimensions")
  }
  s.call <- ds[-MARGIN]
  s.ans <- ds[MARGIN]
  d.call <- d[-MARGIN]
  d.ans <- d[MARGIN]
  dn.call <- dn[-MARGIN]
  dn.ans <- dn[MARGIN]
  d2 <- prod(d.ans)
  newx <- aperm(x, c(s.call, s.ans))
  dim(newx) <- c(prod(d.call), d2)
  ans <- vector("list", d2)
  for (i in seq_len(d2)) {
    ans[[i]] <- array(newx[, i], d.call, dn.call)
  }
  array(ans, d.ans, dn.ans)
}

option_list <- list(
  make_option("--i", default = "AllMet.csv", type = "character", help = "metabolite data file"),
  make_option("--g", default = "SampleInfo.csv", type = "character", help = "sample group file"),
  make_option("--sc", default = "sample_color.txt", type = "character", help = "sample color file"),
  make_option("--config", default = "config.csv", type = "character", help = "config file")
)
opt <- parse_args(OptionParser(option_list = option_list))

configData <- read.csv(opt$config, header = F, stringsAsFactors = F) %>%
  set_colnames(c("arg", "value")) %>%
  column_to_rownames("arg")

configData

isPaired <- configData["isPaired", "value"] == "T"

sampleInfo <- read.csv(opt$g, header = T, stringsAsFactors = F) %>%
  select(c("SampleID", "ClassNote")) %>%
  mutate(ClassNote = as.character(ClassNote)) %>%
  mutate(ClassNote = factor(ClassNote, levels = unique(ClassNote)))

sampleColDf <- read.csv(opt$sc, header = T, stringsAsFactors = F, comment.char = "") %>%
  select(c("ClassNote", "col"))
sampleCols <- sampleColDf %>%
  deframe()

groups <- unique(sampleInfo$ClassNote)

parent <- "./"
fileName <- paste0(parent, "/AllMet_Test.csv")
outFileName <- paste0(parent, "/AllMet_Boxplot_with_Points_and_Significant_labels_post_hoc.pdf")

markerDf <- read_tsv("Markers.txt")

if (nrow(markerDf) == 0) {
  quit(status = 0)
}

pValueData <- read.csv(fileName, header = T, stringsAsFactors = F, comment.char = "") %>%
  arrange(P)

names <- pValueData$Metabolite
data <- read.csv(opt$i, header = T, stringsAsFactors = FALSE) %>%
  select(-c("HMDB", "KEGG", "Class")) %>%
  gather("SampleID", "Value", -Metabolite) %>%
  inner_join(sampleInfo, by = c("SampleID")) %>%
  filter(ClassNote %in% groups)

getPList <- function(name) {
  pData <- data %>% filter(Metabolite == name)
  maxY <- max(pData$Value)
  minY <- min(pData$Value)
  yHeight <- (maxY - minY)
  pValueDf <- pValueData %>% filter(Metabolite == name)
  comparisons <- groups %>%
    as.character() %>%
    combn(2) %>%
    asplit(2)
  compareNum <- comparisons %>%
    length()
  limitY <- c(minY, (maxY + yHeight * 0.1 * compareNum) * 1.025)

  dfMethod <- pValueDf$test.method
  testTb <- if (dfMethod == "anova") {
    test <- aov(Value ~ ClassNote, data = pData)
    testTb <- TukeyHSD(test) %>%
      .$ClassNote %>%
      as.data.frame() %>%
      rownames_to_column("group") %>%
      rowwise() %>%
      do({
        result <- as_tibble(.)
        vs <- result$group %>%
          as.character() %>%
          str_split("-") %>%
          unlist() %>%
          str_trim()
        group1 <- vs[1]
        group2 <- vs[2]
        result %>%
          mutate(group1 = group1, group2 = group2, method = "tukey")
      }) %>%
      ungroup()
    # print(testTb)
    testTb
  }else if (dfMethod == "kruskal.test") {
    testTb <- tryCatch({
      dunnTest(Value ~ ClassNote, data = pData) %>%
        .$res %>%
        as_tibble() %>%
        rename(group = Comparison, `p adj` = P.unadj) %>%
        rowwise() %>%
        do({
          result <- as_tibble(.)
          vs <- result$group %>%
            as.character() %>%
            str_split("-") %>%
            unlist() %>%
            str_trim()
          group1 <- vs[1]
          group2 <- vs[2]
          result %>%
            mutate(group1 = group1, group2 = group2, method = "dunn")
        }) %>%
        ungroup()
    }, error = function(e) {
      tb <- pData$ClassNote %>%
        unique() %>%
        as.character() %>%
        combn(2) %>%
        asplit(2) %>%
        map_dfr(function(groups) {
          tibble(group1 = groups[1], group2 = groups[2], `p adj` = 1, method = "dunn")
        })
      print("=E=")
      print(tb)
      tb
    })

    testTb
  }

  testTb <- testTb %>%
    rename(p = `p adj`) %>%
    mutate(i = 1:n()) %>%
    rowwise() %>%
    do({
      result <- as_tibble(.)
      i <- result$i
      p <- result$p
      y.position <- maxY + yHeight * (i) * 0.1
      p.signif <- if (p < 0.0001) "***" else if (p < 0.01) "**" else if (p < 0.05) "*" else "ns"
      result %>%
        mutate(p.signif = p.signif, y.position = y.position)
    }) %>%
    ungroup() %>%
    filter(p.signif != "ns")

  p <- ggplot(pData) +
    theme_bw(base_size = 8, base_family = "Times") +
    theme(axis.text.x = element_text(size = 8, vjust = 0.5),
          axis.text.y = element_text(size = 8), legend.position = 'none',
          axis.title.y = element_text(size = 11), legend.margin = margin(t = 0.3, b = 0.1, unit = 'cm'),
          legend.text = element_text(size = 6), axis.title.x = element_text(size = 11), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5, size = 8),
          panel.border = element_rect(size = 0.25)
    ) +
    xlab("") +
    ylab("") +
    ggtitle(paste0(pValueDf$Metabolite, "\n", "P=", scientific(pValueDf$P, 2))) +
    stat_boxplot(mapping = aes(x = ClassNote, y = Value), geom = "errorbar", width = 0.5) +
    geom_boxplot(mapping = aes(x = ClassNote, y = Value, fill = ClassNote), outlier.shape = NA) +
    geom_jitter(mapping = aes(x = ClassNote, y = Value, fill = ClassNote), size = 0.75) +
    scale_fill_manual("", values = sampleCols) +
    scale_y_continuous(limits = limitY) +
    stat_pvalue_manual(testTb, label = "p.signif", bracket.size = 0.3, tip.length = 0.015)
  outTestTb <- testTb %>%
    mutate(Metabolite = name)
  list(testTb = outTestTb, p = p)
}

list <- names %>%
  map(getPList)

p <- list %>%
  map(function(l) {
    l$p
  })

testTb <- list %>%
  map_dfr(function(l) {
    l$testTb
  }) %>%
  select(-c("y.position")) %>%
  select(c("Metabolite", "group1", "group2", "p", "p.signif"), everything()) %>%
  select_if(!(names(.) %in% (c("group", "diff", "lwr", "upr", "i", "Z", "P.adj"))))

write_csv(testTb, "post_hoc.csv")

pdf(outFileName, 7.5, 9)

for (i in seq(1, length(p), 9)) {
  plot.new()
  iEnd <- i + 8
  inP <- p[i:iEnd] %>%
    map(function(x) {
      if (is.null(x)) {
        p <- ggplot() + theme_void()
        p
      }else x
    })
  inEnd <- if (length(inP) < 3) {
    length(inP)
  }else 3
  inP[1:inEnd] = inP[1:inEnd] %>%
    map(~.x + theme(plot.margin = margin(t = 0.5, unit = "cm")))
  egg::ggarrange(plots = inP, ncol = 3, newpage = F)
}

dev.off()
